#include<cmath>
#include<iostream>


using namespace std;


class Integrate
{
	double f(double x);
	
public:
	double r;
    double x0,xn;
	double trapezoidal(long n);
	double simpson13(long n);
	double simpson38(long n);
	
};


double Integrate::f(double x)
{
	return sqrt(abs(x*x - r*r));

}


double Integrate::trapezoidal(long n)
{
	/*  nh= xn-x0  */
    x0=0;
	xn=r;
	double h= (xn-x0)/n;

	long i;
	double total=0,tx=x0;
	
	for(i=1;i<n;i++)	
	{
		tx=tx+h;
		total=total+f(tx);
	}	

	double integral= h*(f(x0)+f(xn)+2*total)/2;
	
	return 4*integral/(r*r);

}


double Integrate::simpson13(long n)
{
    x0=0;
    xn=r;
	double h= (xn-x0)/n;
	
	double total4=0,total2=0,tx=x0;
	
	long i;
	for(i=1;i<n;i++)
	{
		tx=tx+h;
		
		if(i%2==0)
		{
			total2=total2+f(tx);
		}
		else
		{
			total4=total4+f(tx);
		}

	}

	double integral= h*(4*total4 + 2*total2 + f(x0) + f(xn))/3;
	
	return 4*integral/(r*r);
	
}


double Integrate::simpson38(long n)
{
    x0=0;
    xn=r;
	double h= (xn-x0)/n;
	
	double total3=0,total2=0,tx=x0;
	
	long i;
	for(i=1;i<n;i++)
	{
		tx=tx+h;
		
		if(i%3==0)
		{
			total2=total2+f(tx);
		}
		else
		{
			total3=total3+f(tx);
		}
	}
		
	double integral= 3*h*(3*total3 + 2*total2 + f(x0) + f(xn))/8;
	
	return 4*integral/(r*r);
}






//  #############        main  starts here       ###################


int main()
{

	cout <<"\n           PI CALCULATION USING NUMERICAL INTEGRATION\n\n";	
	cout.precision(10);
	
	Integrate integrator; 
	
	cout<<"\n Enter R : ";
	cin>>integrator.r;
	
	double ans,lastans;
	long n;
	const double tolerence=0.000001;
	int i;
	
	/*     //test code 
	double t=integrator.trapezoidal(10000);
	double kk=3.5555555555123;
	
	cout<<"\n\nANS : " <<4*t <<"  nothing : " <<kk <<"\n\n";
	*/
	
	n=1;
	i=0;
	ans=integrator.trapezoidal(n);
	do 
	{
		n=n+1;
		lastans=ans;
		ans=integrator.trapezoidal(n);
		i++;
	
	}while(abs(ans-lastans)>tolerence);
	
	cout <<"\nTrapezoidal   : " <<"iterations : " <<i <<"    n= " <<n <<"       ANS : " <<ans;
	
	
	n=2;
	i=0;
	ans=integrator.simpson13(n);	
	do 
	{
		n=n+2;
		lastans=ans;
		ans=integrator.simpson13(n);	
		i++;
	
	}while(abs(ans-lastans)>tolerence);
	
	cout <<"\n\nSimpson's 1/3 : " <<"iterations : " <<i <<"    n= " <<n <<"       ANS : " <<ans;
	
	
	n=3;
	i=0;
	ans=integrator.simpson38(n);	
	do 
	{
		n=n+3;
		lastans=ans;
		ans=integrator.simpson38(n);	
		i++;
	
	}while(abs(ans-lastans)>tolerence);
	
	cout <<"\n\nSimpson's 3/8 : " <<"iterations : " <<i  <<"    n= " <<n <<"       ANS : " <<ans <<"\n\n";
	
	return 0;
}


